Dynamics of Wetting Fronts in Porous Media 
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We propose a new phenomenological approach for describing the dynamics of wetting front prop- 
agation in porous media. Unlike traditional models, the proposed approach is based on dynamic 
nature of the relation between capillary pressure and medium saturation. We choose a modified 
phase-field model of solidification as a particular case of such dynamic relation. We show that in 
the traveling wave regime the results obtained from our approach reproduce those derived from the 
standard model of flow in porous media. In more general case, the proposed approach reveals the 
dependence of front dynamics upon the flow regime. 

PACS: 47.55.Mh, 87.22. As, 47.55.Kf, 68.35. Ja 

The dynamics of fluids in porous media has been a 
subject of numerous theoretical and experimental stud- 
ies because of its importance for engineering and environ- 
mental applications Among the most challenging 
problems in this area is modeling fluid flow through a par- 
tially saturated medium, in particular the propagation 
of wetting (or drying) fronts. This has been addressed 
on both "microscopic" (or pore-scale) and "macroscopic" 
(on scales larger than pore size) levels. The fluid dynam- 
ics in the pore networks has been studied by numerous 
authors both theoretically Q and experimentally |||]]. 
These studies on the microscopic scale provide valuable 
insight into the underlying mechanisms of fluid transport 
in porous media. A microscopic description requires the 
detailed information of the pore structure and pore-size 
distribution. When this information is not available, as 
often happens for large domains, one has to rely upon 
a macroscopic description. This description is useful for 
determining such integral characteristics, as the width 
and propagation velocity of moving wetting fronts, the 
influence of its curvature on the dynamics, etc. 

The subject of the present work is the dynamics of wet- 
ting fronts in porous media when a liquid phase (water) 
displaces air. A straightforward description of this pro- 
cess consists of treating wetting fronts as sharp interfaces 
which separates completely wet and dry regions [§,0. 
However, in many realistic cases the structure of the tran- 
sitional zone (wherein water saturation varies gradually) 
cannot be neglected. One can either describe the dy- 
namics of the liquid and air phases separately, or noting 
that air pressure is close to atmospheric, consider only 
the dynamics of the liquid phase. These two approaches 
generally produce similar results (Ref. [||, p. 213). The 
latter approach seems to be more attractive due to its 
simplicity. 

Within the framework of the chosen approach wa- 



ter flow is described by Darcy's law which, similarly to 
Ohm's law for electric current, stipulates the propor- 
tionality between flux and gradient of a potential (see, 
e.g., [Q). For flow through porous media a water pres- 
sure, averaged over many pore sizes, plays the role of 
the potential. Richards || has empirically generalized 
Darcy's law onto flow in partially saturated porous media 
(PSPM) by letting the proportionality coefficient depend 
on saturation 9 of the medium (0 < 9 < 1). Coupled with 
the mass conservation law, the generalized Darcy's law 
constitutes the existing macroscopic description of the 
fluid dynamics in PSPM. To complete this description, 
one needs to specify a relation between 9 and capillary 
pressure (normalized by the product of fluid density p 
and gravitational acceleration g), ip. Traditionally this 
relation is assumed to be algebraic However, this 

contradicts the numerous experimental evidences reveal- 
ing the dependence of 9 — ip relation upon the conditions 
of the experiment. In particular, this relation exhibits 
hysteresis for wetting and drying of a medium |l],[l0|,[ll] . 

In the present Letter we propose a novel phenomeno- 
logical approach to describe propagation of wetting fronts 
in porous media, which is based on a dynamic 9 — ip re- 
lation. Our description consists of two dynamic equa- 
tions for 9 and ip coupled by nonlinear sources. This im- 
plies that the traditional algebraic 9 — ip correspondence 
is replaced by a nonlocal (integro-differential) relation- 
ship. The proposed model is a new application of the 
well-known "phase-field" approach used to describe so- 
lidification, electro-deposition, and other physical prob- 
lems [^2|-p~7| . We show that, under certain conditions, 
dynamics of the wetting front in our approach reproduce 
that in the traditional approach. The same is true for 
the pressure profiles associated with the wetting fronts 
obtained from both models. We demonstrate that, under 
different conditions, this equivalence breaks down, with 
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our model revealing a dynamic nature of 9 — ip relation. 

The generalized Darcy's law for flux q through PSPM 
has a form q = —K(9)V(%p — X3) , where K is the 
saturation-dependent conductivity of the medium and 
X3 is the vertical coordinate (positive downward) that 
stands for the gravitational component of pressure (nor- 
malized by pg). Here ip is measured relative to the am- 
bient atmospheric pressure, so that ip < for partially 
saturated and ip > for fully saturated medium. Com- 
bined with the mass conservation law, d9/dt — — V • q, 
this leads to the equation introduced by Richards pt] 
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Over the years a number of algebraic — ip and K(9) rela- 
tions have been proposed empirically (see, e.g. Jl8|-p0t). 
Richards equation (RE) ([j]) combined with such a rela- 
tion constitutes the traditional approach to describe flow 
in PSPM. Since algebraic relations are not supported by 
experimental data, in what follows we propose a new ap- 
proach with a dynamic 9 — ip relation. 

To describe flow in PSPM, we modify the phase-field 
model (PFM) from Q 
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where D = K s /S is the diffusion coefficient, K s is the 
conductivity of a fully saturated medium, S is the spe- 
cific storage (measure of compressibility of the fluid and 
medium), r is a characteristic time-scale of the satura- 
tion dynamics, and W is the width of a moving wetting 
front. The model parameter A will be determined as a 
function of macroscopic parameters K s , W and r . The 
constant ipf is the normalized capillary pressure along 
the moving front in the sharp-interface limit. We will 
demonstrate below that ipf coincides with the pressure 
in the dry medium far ahead of the wetting front. Since 
the width of the capillary zone, associated with the local- 
ized wetting fronts, is much smaller than the typical scale 
of the pressure variation, the following holds W 2 /t <C D. 
We have added the last term in (§) to PFM Q, to ac- 
count for the gravitational force. 

Our model is phenomenological in the sense that 
presently we do not provide a rigorous physical motiva- 
tion for the nonlinear source term on the right-hand-side 
of (^). Nevertheless, the model captures the main fea- 
tures of the wetting fronts propagation in porous media. 
In particular, numerous experiments |^],^2| have shown 
that, under certain conditions, the wetting fronts remain 
localized and propagate in a self-similar manner. The 
medium is fully saturated (9=1) behind the front re- 
gion and completely dry (9 = 0) ahead of the front. The 



cubic polynomial in (g) provides for such a structure of 
the wetting front (similar to the PFM for solidification). 
The fact that the liquid moves in the direction opposite 
to the pressure gradient is accounted for by the term pro- 
portional to (ip — ipf) in (||), since its presence makes the 
depths of the two minima of the corresponding potential 
energy different. 

Eqs. (§)-(||) have to be supplemented by a constraint 
to ensure conservation of mass. The specific storage is re- 
lated to the compressibility of fluid and porous medium 
as S = pguj(P f - p s + P ) (Ref. @, p. 108). Here 
u> is the medium porosity (fraction of pore volume in 
the total volume v of the medium) and (3f,j3 m , and (3 P 
are the compressibility coefficients of fluid, solid grains, 
and pores, respectively. The total mass of the fluid is 
given by M = f y uip9dx. For incompressible fluids 
and media (S = 0), since D = K s /S, Eq. (||) reads 
89/dt = K s (W 2 ip — 89/8x3) , which conserves mass. For 
compressible fluids and media dM/dt = Q, where Q is 
the total mass flux through the medium, provides the 
global constraint on the system (||)-(|]). 

Since S is the specific storage of a medium, the effect 
of compressibility of fluids and media on flow dynam- 
ics is characterized by the dimensionless parameter SL, 
where L is a typical domain size. For many practical 
applications, such as flow of water through low poros- 
ity rocks, S - 1CT 5 - lCT 7 ™" 1 and L - 1CT 1 - 10 3 m. 
Thus SL -C 1, and the fluid and medium are virtually 
incompressible. 

We consider propagation of one-dimensional (ID) wet- 
ting fronts in a slightly compressible medium with SL <C 
1. The dynamics of these fronts is described by ID Eqs. 
(|)-(|), and the mass conservation constraint is satisfied 
automatically. Two different situations exist: horizontal 
(gravity-free) and vertical front propagation. To main- 
tain propagation of a self-similar front, we choose the con- 
stant flux condition d(ip — xs)/dxi = —q/K s at Xi = 0, 
and no-flux condition for capillary force dip/dxt = at 
Xi = L. Here i = 1 or i = 3 for horizontal or vertical 
flow, respectively, and q is an external velocity flux. At 
the boundary behind the front, Xi = 0, the medium is 
fully saturated (9 = 1), and at the boundary ahead of 
the front, Xi = L, the medium is dry (9 = 0). 

Let us apply a traveling wave ansatz z = X{ — Vt, for 
a front moving with velocity V, to ID Eqs. (§)-(|). Re- 
lations W ~ Vt and W 2 /t <C D give rise to a small 
parameter Pe <C 1, where Pe = VW/D is the Peclet 
number. The perturbation analysis around Pe = 0, sim- 
ilar to that performed in Jl^ ], gives a front moving with 
velocity V = q, and saturation and capillary pressure 
profiles 
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Here J is an absolute value of the capillary pressure gra- 
dient at the boundary. It is given by J — q/K s for a hor- 
izontally propagating front, and by J = (q/K s — 1) for 
a vertically propagating front. A vertical wetting front 
does not exist when q < K s , i.e. when the external flux is 
not large enough to compensate the effect of the gravity 
force, and thus develop a saturated zone. 

The saturation profile (0) is evaluated in the zeroth 
order in Pe, since for the localized wetting fronts satura- 
tion varies in the narrow region W <C V Dt. Hence the 
higher order corrections to 9 influence the nonlocal pres- 
sure profile (0) only in the second order. The solvability 
condition for the equation for 9 in the first order (similar 
to p|) yields r/A « -0.313W 2 J/q . 
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FIG. 1. Normalized capillary pressure (ip — ipf)K s /qW 
(a) and saturation (b) profiles in traveling wave coordinate. 
In (b) normalized capillary pressure in a dry medium is taken 
ip f K a /qW = -10. 

We now compare our approach with the traditional 
(Richards) approach. One of the most widely-used ex- 
amples of algebraic functions 8(ip) and K(9), complemen- 
tary to RE (pi), has been proposed by Gardner [^0|Jl^1 



0WO 



K{9) = K s 9 
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where a is the pore-size distribution parameter. The first 
equation in (^j) is valid for ip < (when the medium is 
partially saturated), while for ip > (for a fully satu- 
rated medium) 9 = 1. The same boundary conditions as 
before are used. Substituting (^) into ID Eq. (|l|) and 
applying the traveling wave ansatz gives, after a series of 
transformations, the solution for the pressure profile 
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where ipL is a large negative number corresponding to 
the pressure in the dry medium far ahead of the front. 

Comparing (Q) with (||) we find that the pressure pro- 
files obtained from both models coincide, provided the re- 
lations between parameters \/2WJ — a^ 1 and ipf = ipL 
hold. Figure [j](a) shows the resulting pressure profile 
from both models. Substituting ^ into the first rela- 
tion in ([}]) we obtain the saturation profile for Richards 
model. The result is compared with the saturation profile 
([|) in the Figure 0(b). Although the saturation profiles 
differ within the capillary zone, this difference occurs on 
the very small scale W . An unphysical sharp kink that 
appears in the Richards model for 9(z) is absent in our 
model. Comparing Figure 0(b) with the experimental 
data presented on Figures 2-5 of Ref. Q shows that the 
saturation profile obtained by our model agrees with the 
experiments. 

Though in the traveling wave regime our model repro- 
duces the results obtained from the Richards model, it 
is also capable of capturing the dynamic 9 — ip relation 
in different regimes, without adjusting the parameters. 
This is not the case for the Richards model, where one 
needs to adjust the parameter a in (^|) to different ex- 
perimental regimes. To demonstrate the dynamic nature 
of the 9 — ip relation in our model, we have performed 
the numerical simulations of ID Eqs. (||)-(||) under sev- 
eral boundary conditions. We solved system (§)-(§]) us- 
ing finite differences. Note that though in our simula- 
tions SL <C 1, we keep the term Sdip/dt in (g). Despite 
being small, this term provides relaxational feature to 



numerical solution of (a 
cal algorithm. Figure p2 



which stabilizes the numeri- 
shows 9(tp) obtained from these 
simulations for horizontal flow. The solid line in the fig- 
ure corresponds to the traveling wave regime, resulting 
from the previously described boundary conditions. The 
dashed line represents 9(ip) at time t — 1000 correspond- 
ing to the wetting with the constant pressure at x\ = 
and the same no-flux condition at X\ = L . 

Note that in the limit of the narrow capillary zone, 
W — ► , our model describes the dynamics of sharp 
interface which separates completely wet (0 = 1) and 
completely dry [9 — 0) regions. This is in analogy 
with the PFM of solidification that reduces to the free- 
boundary problem Jl2"|| . It follows from (||) that in this 
limit ip(0) — ipf and ip z (0) = —J. Moreover, (||) re- 
duces to a diffusion equation which is commonly used to 
describe flow in saturated media. 

In conclusion, we have proposed a novel approach to 
describe dynamics of wetting fronts in porous media. Un- 
like the traditional approach, our phenomenological ap- 
proach reflects the dynamic nature of the relation be- 
tween capillary pressure and saturation of a medium. 
We have found that this relation varies with the flow 
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regimes, which is supported by experimental data. We 
have demonstrated that, in the traveling wave regime, 
the proposed model reproduces the results obtained from 
the standard approach. We plan to extend our novel ap- 
proach to two- and three-dimensional media incorporat- 
ing in the description the front curvature and anisotropy 
of the medium. We expect to develop experimental sup- 
port for the proposed approach by measuring quantita- 
tive features of the wetting fronts in porous media, such 
as the dependence of the front width on external flux. 
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FIG. 2. Dynamic 6 — tp relation for both traveling wave and 
non-traveling wave regime. Parameters of the simulations are 
D = 1,W = 0.1,S = 0.005 , r = 1 . For the traveling wave 
regime q — K s . Value of A is calculated according to the 
selection result r/X « — 0.313W 2 J/q . System size L = 10 
and number of grid-points is 201. Time step dt = 0.001 . 
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